1. Coping with Mortality

Learning Objectives

  • Understand options for modeling background mortality based on life tables and parametric mortality models.
  • Construct a cause-deleted life table to model cause-specific and non-cause-specific death.

Approaches to Modeling Mortality

  • Standard approach is to draw from life table data:
age lx dx mx qx
0 100000 553.22 0.00556 0.00553
1 99447 38.18 0.00038 0.00038
2 99409 23.16 0.00023 0.00023
  • lx is the total number living in a hypothetical cohort (of 100,000).
  • dx is the number of deaths in the age interval.

Approaches to Modeling Mortality

  • Standard approach is to draw from life table data:
age lx dx mx qx
0 100000 553.22 0.00556 0.00553
1 99447 38.18 0.00038 0.00038
2 99409 23.16 0.00023 0.00023
  • mx is the death rate within the interval.
  • qx is the probability of death within an interval.

Life Table Data

  • Death probabilities are useful for discrete time Markov models.
age lx dx mx qx
0 100000 553.22 0.00556 0.00553
1 99447 38.18 0.00038 0.00038
2 99409 23.16 0.00023 0.00023

Life Table Data

  • Death rates are useful for discrete event simulation (DES) models.
age lx dx mx qx
0 100000 553.22 0.00556 0.00553
1 99447 38.18 0.00038 0.00038
2 99409 23.16 0.00023 0.00023

Life Table-Based Mortality

  • A common approach is to use a life table-based lookup table in the model:
tpDn_lookup <-
    c("(34,44]" = 0.0017,  
      "(44,54]" = 0.0044,
      "(54,64]" = 0.0138,
      "(64,74]" = 0.0379,
      "(74,84]" = 0.0912,
      "(84,100]" = 0.1958)
  • Highlighted row: annual probability of death for 35-44 year old.

Life Table-Based Mortality

  • What if age bins are coarse (e.g., 0, 5, 10,… vs. 0,1,2,3,…)?
    • Assume death probability/rate is constant within age bins (and cycles)?
    • In a DES model, we need to sample a time to death.
  • Can sample from survival function in underlying life table.
  • Ages in life table often topcoded (e.g., 85, 100, etc.)

Alternative: Parametric Mortality Model

  • Fit a parametric mortality model to the life table data.
  • Reduces dimensionality to a few parameters.
  • Use these parameters to sample time of death, obtain mortality rate/probability at any arbitrary time, etc.

Obtaining Life Table Data

The R packages MortalityLaws and demography facilitate life table data aquisition and modeling.

Common sources:

Processing Life Table Data

  • Code below extracts U.S. life table data from the human mortality database (HMD).

Processing Life Table Data

  • Code below extracts U.S. life table data from the human mortality database (HMD).
hmd_usa <-
    demography::hmd.mx("USA",
                       username = "<<your user name>>",
                       password = "<<your password>>", "USA")

Processing Life Table Data

  • Code below processes HMD data and constructs a (2019) life table from it.
  • Radix (cohort size) is 100,000 individuals.

Processing Life Table Data

  • Code below processes HMD data and constructs a (2019) life table from it.
  • Radix (cohort size) is 100,000 individuals.
mortality_year = 2019 
radix = 100000

lt = 
    hmd_usa %>% 
    demography::lifetable(.,series = "total", years = mortality_year) %>% 
    as_tibble() %>% 
    rename(age = x) %>% 
    mutate_at(vars(lx,dx), function(x) x * radix) 

US Life Table (2019)

Highlighted column: Remaining life expectancy at exact age x.

age lx dx ex qx mx
0 100000 553.221 78.963 0.00553 0.00556
1 99447 38.180 78.402 0.00038 0.00038
2 99409 23.160 77.432 0.00023 0.00023
3 99385 17.689 76.450 0.00018 0.00018

Alive-Dead Model

Alive-Dead Model

  • As a refresher, let’s construct a very simple alive-dead model using life table death probability (qx) values.
G Alive Alive Alive->Alive 1.0-qx(t) Dead Dead Alive->Dead qx(t) Dead->Dead 1.0

Parameterize the model

params = 
  list(
    t_names = c("lifetable"),         # Strategy names. 
    n_treatments = 1,                 # Number of treatments
    s_names  = c("Alive", "Dead"),    # State names
    n_states = 2,                     # Number of states
    n_cohort = 100000,                # Cohort size
    initial_age = 0,                  # Cohort starting age    
    n_cycles = 110,                   # Number of cycles in model.  
    cycle = 1,                        # Cycle length
    life_table = lt                   # Processed HMD life table data (2019, USA)
  )

Parameterize the model

params = 
  list(
    t_names = c("lifetable"),         # Strategy names. 
    n_treatments = 1,                 # Number of treatments
    s_names  = c("Alive", "Dead"),    # State names
    n_states = 2,                     # Number of states
    n_cohort = 100000,                # Cohort size
    initial_age = 0,                  # Cohort starting age    
    n_cycles = 110,                   # Number of cycles in model.  
    cycle = 1,                        # Cycle length
    life_table = lt                   # Processed HMD life table data (2019, USA)
  )
  • No specific strategies or policies under consideration; we’re just trying to replicate mortality using life table data for a hypothetical cohort of 100,000 infants.

Parameterize the model

params = 
    list(
        t_names = c("lifetable"),         # Strategy names. 
        n_treatments = 1,                 # Number of treatments
        s_names  = c("Alive", "Dead"),    # State names
        n_states = 2,                     # Number of states
        n_cohort = 100000,                # Cohort size
        initial_age = 0,                  # Cohort starting age    
        n_cycles = 110,                   # Number of cycles in model.  
        cycle = 1,                        # Cycle length
        life_table = lt                   # Processed HMD life table data (2019, USA)
    )
  • Discrete time will be in annual cycles.

Parameterize the model

params = 
    list(
        t_names = c("lifetable"),         # Strategy names. 
        n_treatments = 1,                 # Number of treatments
        s_names  = c("Alive", "Dead"),    # State names
        n_states = 2,                     # Number of states
        n_cohort = 100000,                # Cohort size
        initial_age = 0,                  # Cohort starting age    
        n_cycles = 110,                   # Number of cycles in model.  
        cycle = 1,                        # Cycle length
        life_table = lt                   # Processed HMD life table data (2019, USA)
    )
  • The lifetable “strategy” will draw on the HMD life table data processed earlier.

Age-Dependent Transition Matrix

  • We’ll next construct a transition probability matrix using the age-specific probability of death field (qx) from the life table data.

  • To do so, we need to define a function that, for a given age and cycle length, calculates the probability of death and places this within a transition matrix.

Age-Dependent Transition Matrix

Life Table:

lt[1,c("age","qx")]
# A tibble: 1 × 2
    age      qx
  <dbl>   <dbl>
1     0 0.00553
lt[90,c("age","qx")]
# A tibble: 1 × 2
    age    qx
  <dbl> <dbl>
1    89 0.123

Transition Matrix:

params$mP[[1]][,,"lifetable"]
       to
from    Alive   Dead     
  Alive 0.99447 0.0055322
  Dead  0       1        
params$mP[[90]][,,"lifetable"]
       to
from    Alive   Dead   
  Alive 0.87677 0.12323
  Dead  0       1      

Construct the Markov Trace

sim_cohort <- function(params) {
 params$t_names %>% map(~({ 
    tr_ <- t(c("alive" = params$n_cohort, "dead" = 0))
    
    res <- do.call(rbind,lapply(params$mP, function(tp) {
        tr_ <<- tr_ %*% matrix(unlist(tp[,,.x]),nrow=params$n_state)
    }))
    res <- rbind(c(params$n_cohort,0),res) 
    dimnames(res) <- list(paste0(c(0:params$n_cycles)), params$s_names)
    res
  })) %>% 
    set_names(params$t_names)
}

Construct the Markov Trace

  • Output: named list object with full Markov trace for each state.

Construct the Markov Trace

  • Output: named list object with full Markov trace for each strategy.
trace <- 
  sim_cohort(params)
trace[["lifetable"]][1:10,]
   Alive   Dead
0 100000   0.00
1  99447 553.22
2  99409 591.40
3  99385 614.56
4  99368 632.25
5  99354 646.46
6  99340 659.67
7  99328 671.69
8  99317 682.52
9  99307 693.14

Checkpoint: Life Expectancy at Birth

Life Expectancy

  • We have a model with 110 annual cycles for a cohort of 100,000 0 year olds.
  • According to the underlying life table data, life expectancy at birth is 78.96315 years.
  • Let’s verify we can replicate this with our Markov model.

Life Expectancy

  • We have a model with 110 annual cycles for a cohort of 100,000 0 year olds.
  • According to the underlying life table data, life expectancy at birth is 78.96315 years.
  • Let’s verify we can replicate this with our Markov model.
age qx ex
0 0.00553 78.963

Life Expectancy

  1. Define the payoff for life expectancy (1 if alive, 0 otherwise)
# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

Life Expectancy

  1. Calculate life years in each cycle.1
# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

Life Expectancy

  1. Define a function for cycle adjustment.1
# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

Life Expectancy

  1. Sum of values (with cycle adjustment applied) is total life expectancy.
# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

Life Expectancy

Age Life Expectancy (Life Table data) Life Expectancy (Markov-Life Table)
0 78.963 78.924
  • We can repeat this exercise for different cohort starting ages and verify things line up.

Life Expectancy

Mortality Modeling

Mortality Modeling

  • Depending on context and modeling environment, you may need to model mortality.
  • The MortalityLaws package has a number of mortality models we can draw from:

Mortality Modeling

Name Model Code
Infant mortality
Opperman \(\mu[x] = A/\sqrt{x} - B + C\sqrt{x}\) opperman
Weibull \(\mu[x] = 1/\sigma * (x/M)^{(M/\sigma - 1)}\) weibull
Accident hump
Inverse-Gompertz \(\mu[x] = [1-\exp((M-x)/\sigma)]/[\exp((M-x)/\sigma)-1]\) invgompertz
Inverse-Weibull \(\mu[x] = 1/\sigma * (x/M)^{[-M/\sigma - 1]} / [\exp((x/M)^{(-M/\sigma)}) - 1]\) invweibull

Name Model Code
Adult mortality
Gompertz \(\mu[x] = A \exp[Bx]\) gompertz
Gompertz \(\mu[x] = 1/\sigma \exp[(x-M)/\sigma]\) gompertz0
Makeham \(\mu[x] = A \exp[Bx] + C\) makeham
Makeham \(\mu[x] = 1/\sigma \exp[(x-M)/\sigma)] + C\) makeham0
Perks \(\mu[x] = [A + BC^x] / [BC^{-x} + 1 + DC^x]\) perks
Strehler-Mildvan \(\mu[x] = K \exp[-V_0 (1-Bx)/D]\) strehler_mildvan
Adult and/or old-age mortality
Van der Maen \(\mu[x] = A+Bx+Cx^2 + I/[N-x]\) vandermaen
Beard \(\mu[x] = \exp(Bx)/[1+KA \exp(Bx)]\) beard
Beard-Makeham \(\mu[x] = A \exp(Bx)/[1 + KA \exp(Bx)]\) beard_makeham
Gamma-Gompertz \(\mu[x] = A \exp(Bx)/(1+AG/B[\exp(Bx)-1])\) ggompertz

Name Model Code
Old-age mortality
Van der Maen \(\mu[x] = A + Bx + I/[N-x]\) vandermaen2
Quadratic $= A + Bx + Cx^2% quadratic
Kannisto \(\mu[x] = A \exp(Bx)/[1+A\exp(Bx)]\) kannisto
Kannisto-Makeham \(\mu[x] = A \exp(Bx)/[1+A\exp(Bx)]+C\) kannisto_makeham

Name Model Code
Full age range
Thiele \(\mu[a] = A e^{-Bx} + C e^{-\frac{1}{2} D(x-E)^2}+F e^{Gx}\) thiele
Wittstein \(q[x] = (1/B) A^{-[(Bx)^N]} + A^{-[(M-x)^N]}\) wittstein
Heligman-Pollard \(q[x]/p[x] = A^{[(x + B)^C]} + D e^{-E \log(x/F)^2} + G H^x\) HP
Heligman-Pollard \(q[x] = A^{[(x + B)^C]} + D e^{-E \log(x/F)^2} + GH^x / [1 + GH^x]\) HP2
Rogers-Planck \(q[x] = A_0 + A_1 e^{-Ax} + A_2 e^{B(x-u)-e^{-C(x-u)}}+A_3 e^{Dx}\) rogersplanck
Martinelle \(\mu[x] = [A e^{Bx} + C]/[1+D e^{Bx}]+K e^{Bx}\) martinelle
Carriere \(l[x] = P1 l[x](weibull) + P2 l[x](invweibull) + P3 l[x](gompertz)\) carriere1
Kostaki \(q[x]/p[x] = A^{[(x+B)^C]} + D \exp[-(E_i \log(x/F_))^2] + G H^x\) kostaki

Mortality Modeling

  • Another nice feature of the package is that each mortality model type has its own defined function to calculate the mortality rate (or probability) given an age and model coefficients.

Mortality Modeling

Generally speaking, we need three inputs:

  • age: Ages for lifetable
  • dx: The number of deaths between exact ages x and x+1.
  • lx: Number of survivors to exact age x.

Gompertz Model

  • Because we want to stay general (i.e., model all over the age spectrum), our first attempt will be a Gompertz model.
  • Gompertz reduces mortality to two parameters.

Gompertz Model

ages     <- lt$age[lt$age<100]
deaths   <- lt$dx[lt$age<100]
exposure <- lt$lx[lt$age<100]

gom_fit <- MortalityLaw(
    x  = ages,      # vector with ages
    Dx  = deaths,   # vector with death counts
    Ex  = exposure, # vector containing exposures
    law = "gompertz",
    opt.method = "LF2")

Gompertz Model

# Fitted coefficients
gom_fit$coefficients
         A          B 
0.00010772 0.07534542 

Gompertz Model

# Fitted coefficients
gom_fit$coefficients
         A          B 
0.00010772 0.07534542 
# Hazard rate of death for a 40 year old based on Gompertz model
# mx = A exp(Bx)
gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )
        A 
0.0021938 

Gompertz Model

# Fitted coefficients
gom_fit$coefficients
         A          B 
0.00010772 0.07534542 
# Hazard rate of death for a 40 year old based on Gompertz model
# mx = A exp(Bx)
gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )
        A 
0.0021938 
# Can simply use the supplied function to calcualte 
MortalityLaws::gompertz(x = 40, par = gom_fit$coefficients)
$hx
[1] 0.0021938

$par
         A          B 
0.00010772 0.07534542 

$Sx
[1] 0.97269

Gompertz Model

# Calculate death hazard rate at age 40 from fitted Gompertz model. 
mx = gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

Gompertz Model

# Calculate death hazard rate at age 40 from fitted Gompertz model. 
mx = gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

# Convert mortality rate to probability
-log(1-mx)
        A 
0.0021962 

Gompertz Model

# Calculate death hazard rate at age 40 from fitted Gompertz model. 
mx = gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

# Convert mortality rate to probability
-log(1-mx)
        A 
0.0021962 
# Compare with life table probability (qx) at age 40:
lt[41,c("age","qx")] 
# A tibble: 1 × 2
    age      qx
  <dbl>   <dbl>
1    40 0.00208

Gompertz Model

  • Gompertz doesn’t fit all age ranges well—particularly young ages.
  • (This is a well-known fact in demography.)
  • If we were to focus only on adults, however, this would be a nice way to go.

Gompertz Model

plot(gom_fit)

Gompertz Model

  • If we were to focus only on adults, however, this would be a nice way to go.

Gompertz Model: 40 Year Old Cohort

ages     <- lt$age[lt$age>=40 & lt$age<100]
deaths   <- lt$dx[lt$age>=40 & lt$age<100]
exposure <- lt$lx[lt$age>=40 & lt$age<100]

gom_fit40 <- MortalityLaw(
                x  = ages,      # vector with ages
                Dx  = deaths,   # vector with death counts
                Ex  = exposure, # vector containing exposures
                law = "gompertz",
                opt.method = "LF2")

Gompertz Model: 40 Year Old Cohort

plot(gom_fit40)

Heligman-Pollard

  • While you could draw from any of the aforementioned mortality models, the Heligman-Pollard model tends to fit the entire age distribution reasonably well.

  • Easy to use: just use HP in lieu of gompertz!

Heligman-Pollard

  • While you could draw from any of the aforementioned mortality models, the Heligman-Pollard model tends to fit the entire age distribution reasonably well.

  • Easy to use: just use HP in lieu of gompertz!

ages     <- lt$age[lt$age<100]
deaths   <- lt$dx[lt$age<100]
exposure <- lt$lx[lt$age<100]

hp_fit <- MortalityLaw(
                x  = ages,      # vector with ages
                Dx  = deaths,   # vector with death counts
                Ex  = exposure, # vector containing exposures
                law = "HP",
                opt.method = "LF2")

Heligman-Pollard

plot(hp_fit)

Heligman-Pollard

Mortality rate for a 40 year old:

HP(x = 40, par = hp_fit$coefficients)
$hx
[1] 0.0020409

$par
           A            B            C            D            E           F_ 
 0.000386408  0.021271772  0.107422357  0.001025285  2.871760263 33.201213579 
           G            H 
 0.000022277  1.102611045 
x   = 40
mu1 = A^((x + B)^C) + G * H^x)
mu2 = D * exp(-E * (log(x/F_))^2))
eta = ifelse(x == 0, mu1, mu1 + mu2)
hx  = eta/(1 + eta)

Death Hazard at Age 40: Gompertz vs. HP

gompertz(x = 40, par = gom_fit$coefficients)$hx
[1] 0.0021938
HP(x = 40, par = hp_fit$coefficients)$hx
[1] 0.0020409

Alive-Dead (Modeled Mortality)

We’ll Approach This From Two Angles

  1. Discrete Event Simulation: Sample a background mortality time.
  2. Discrete Time Markov: Use parameterized mortality in Alive-Dead model.

1. Sampling Death Dates (DES)

  • Difficult to sample death times from life tables due to coarseness of data (e.g., five-year age bins).
  • Solution:
    1. Fit a parametric mortality model (e.g., Heligman-Pollard)
    2. Construct a survival function from modeled hazards.
    3. Sample a death time from the survival function.

1. Sampling Death Rates (DES)

  • Intuition is easier if we start with a defined quantile (e.g., 25th percentile of survival)

1. Sampling Death Rates (DES)

  • Intuition is easier if we start with a defined quantile (e.g., 25th percentile of survival)
u = 0.25

1. Sampling Death Rates (DES)

  • We can feed the quantile (0.25) and coefficients into an inverse quantile function.
  • Straightforward to do for Gompertz model
u = 0.25
qgompertz(p = u, shape = gom_fit$coefficients["B"], rate = gom_fit$coefficients["A"])
[1] 70.467

1. Sampling Death Rates (DES)

  • Let’s verify this gets us close to what we see in the life table data.
u = 0.25
qgompertz(p = u, shape = gom_fit$coefficients["B"], rate = gom_fit$coefficients["A"])
[1] 70.467
age mx Sx
70 0.01827 0.76590
71 0.01995 0.75077
72 0.02162 0.73472

1. Sampling Death Rates (DES)

  • There is no inverse quantile function for Heligman-Pollard, but we can create one!

1. Sampling Death Rates (DES)

  • There is no inverse quantile function for Heligman-Pollard, but we can create one!
  • Approxfun returns a function that provides a value of y for a given value of x.
qHP <- 
  approxfun(y = 0:115,
            x = 1 - exp(-cumsum(HP(x = 0:115, par = hp_fit$coefficients)$hx)))

1. Sampling Death Rates (DES)

  • For approxfun():
    • x is quantile in CDF for survival, i.e., 1 - Survival1
    • y is the age.
  • So qHP(0.5) would return median survival age in the modeled data.

1. Sampling Death Rates (DES)

qHP <- 
  approxfun(y = 0:115,
            x = 1 - exp(-cumsum(HP(x = 0:115, par = hp_fit$coefficients)$hx)))
qHP(0.25)
[1] 71.069

Survival in life table:

age mx Sx
70 0.01827 0.76590
71 0.01995 0.75077
72 0.02162 0.73472

1. Sampling Death Rates (DES)

  • We now have the tools we need to sample death times that closely match underlying mortality data!
  • For a given simulated patient, sample a uniform (0,1) and feed this number into the inverse quantile function.

1. Sampling Death Rates (DES)

  • We now have the tools we need to sample death times that closely match underlying mortality data!
  • For a given simulated patient, sample a uniform (0,1) and feed this number into the inverse quantile function.
# Sample death times for 10 individuals
u <- runif(10, min = 0, max = 1)
qHP <- 
  approxfun(y = 0:115,
  x = 1 - exp(-cumsum(HP(x = 0:115, par = hp_fit$coefficients)$hx)))
qHP(u)
 [1] 84.613 77.082 87.010 78.662 86.211 92.977 68.957 89.789 64.050 86.566

2. Modeled Mortality (Markov)

  • Let’s now incorporate modeled mortality parameters into our Alive-Dead model.
  • Rather than draw a probability of death from the life table data, we’ll calculate the probability of death at a given age using the model coefficients.

Parameterize the model

params = 
  list(
    t_names = c("lifetable",          # Strategy names
                "heligman-pollard",
                "gompertz"),          
    n_treatments = 3,                 # Number of treatments
    s_names  = c("Alive", "Dead"),    # State names
    n_states = 2,                     # Number of states
    n_cohort = 100000,                # Cohort size
    initial_age = 40,                 # Cohort starting age    
    n_cycles = 110,                   # Number of cycles in model.  
    cycle = 1,                        # Cycle length
    life_table = lt,                  # Processed HMD life table data (2019, USA)
    gom_fit = gom_fit,                # Fitted Gompertz model object
    hp_fit = hp_fit                   # Fitted Heligman-Pollard model object
  )

Age-Dependent Transition Matrices

fn_mPt <- function(t, params) {
  with(params, {
    h = cycle # Time step
    lapply(t, function(tt){
       # Life table method
        current_age_lt = pmin(initial_age  + (tt)*h - 1,max(lt$age))
        p_death <- lt[lt$age==current_age_lt,"qx"]
        # parametric mortality methods
        current_age = initial_age  + (tt)*h - 1
        r_death_gompertz <- gompertz(current_age,params$gom_fit$coefficients)$hx
        p_death_gompertz <- 1 - exp(-r_death_gompertz)
        r_death_hp <- HP(current_age, params$hp_fit$coefficients)$hx
        p_death_hp <- 1 - exp(-r_death_hp)
        
        array(data = c(1 - p_death, 0, 
                 p_death,1,
                 
                 1 - p_death_hp, 0, 
                 p_death_hp,1,
                 
                1 - p_death_gompertz, 0, 
                 p_death_gompertz,1),
              
              dim = c(n_states, n_states, n_treatments),
              dimnames = list(from = s_names,
                              to = s_names,
                              t_names)) 
    })
  })
}

Construct a Markov Trace

# Add transition matrices into params.
params$mP <- 
  fn_mPt(1:params$n_cycles, params)

# Construct the Markov trace. 
trace <- 
  sim_cohort(params)

Construct a Markov Trace

trace[["lifetable"]][1:10,]
   Alive    Dead
0 100000    0.00
1  99792  207.68
2  99584  416.43
3  99365  634.88
4  99138  861.67
5  98899 1101.19
6  98639 1361.45
7  98358 1641.68
8  98054 1945.62
9  97728 2271.99
trace[["gompertz"]][1:10,]
   Alive    Dead
0 100000    0.00
1  99781  219.14
2  99545  454.89
3  99292  708.46
4  99019  981.15
5  98726 1274.35
6  98410 1589.51
7  98072 1928.21
8  97708 2292.12
9  97317 2682.98
trace[["heligman-pollard"]][1:10,]
   Alive    Dead
0 100000    0.00
1  99796  203.88
2  99584  416.04
3  99363  637.42
4  99131  869.10
5  98888 1112.31
6  98632 1368.43
7  98361 1639.00
8  98074 1925.75
9  97769 2230.56

Calculate Life Expectancy

# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = lapply(trace, function(tr) (tr * (1 / params$n_cohort)) %*% payoff_life_exp)

# Function for cycle correction (alternative Simpson's method)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = lapply(life_exp_cycle,cycle_adj)

Life Expectancy: Alive-Dead Model

total_life_exp
$lifetable
[1] 40.959

$`heligman-pollard`
[1] 41.008

$gompertz
[1] 41.29

Cause-Deleted Mortality

Motivation

  • Often, background mortality is drawn directly from life table data.
  • Models also frequently include cause-specific death as a separate tracked event or health state.

Motivation

  • This approach “double-counts” death because cause-specific death is also reflected in the life table estimates!
  • OK if cause-specific death is a small contributor to overall death rates.
  • But what if modeled disease is a frequent cause of death (e.g., cardiovascular disease, cancer)?

Solution

  • We can construct an alternative life table that “nets out” cause-specific deaths.
  • Not foolproof: need to make an assumption that we can parse out deaths as indepdendent contributors to overall mortality.

CVD Natural History Model

CVD Natural History Model

  • Next set of slides will walk through the process of constructing a natural history model for cardiovascular disease.
  • We’ll construct a model with non-CVD background mortality, and CVD mortality.
  • First step is constructing a cause-deleted life table.
    • Needed because at older ages, CVD is a substantial contributor (~50%) of overall mortality.

CVD-Deleted Life Table

Constructing the cause-deleted life table requires two necessary inputs:

  1. Overall mortality by age.
  2. Cause-specific mortality by age (ideally as a % of all deaths at a given age).
  • These are obtainable for many countries from the Global Burden of Disease and Human Cause of Death Database websites.

CVD-Deleted Life Table

Screenshot of Global Burden of Disease Data

We extract the percentage of overall deaths that are from CVD by age bin:

ihme_cvd <- 
  tibble::tribble(
        ~age_name,        ~val,
               1, 0.038771524,
               5, 0.038546046,
              10, 0.044403585,
              15, 0.033781126,
              20, 0.035856165,
              25, 0.053077797,
              30, 0.086001439,
              35, 0.130326551,
              40, 0.184310334,
              45,  0.21839762,
              50, 0.243705394,
              55, 0.256334637,
              60,  0.26828001,
              65, 0.272698709,
              70,  0.28529754,
              75, 0.310642009,
               0, 0.016750489,
              80, 0.353518012,
              85, 0.399856716,
              90, 0.447817792,
              95, 0.495305502
        ) %>% 
    mutate(age_ihme = cut(age_name,unique(c(0,1,seq(0,95,5),105)),right=FALSE))  %>% 
    select(age_ihme,  pct_cvd = val) 

Merge these percentages into the underlying life table data and use them to calculate the total number of CVD deaths by age:

age age_ihme pct_cvd dx dx_i
0 [0,1) 0.01675 553.221 9
10 [10,15) 0.04440 11.518 1
25 [25,30) 0.05308 103.706 6
50 [50,55) 0.24371 365.721 89
75 [75,80) 0.31064 1997.766 621
98 [95,105) 0.49531 1262.593 625

We next calculate the cause-deleted probability of death (qd) and death rate (md), as well as the cause-specific probability of death (qi) and death rate (mi):

Steps

  1. Cause-specific probability of death (qi): \(q \cdot dx_i / dx\)
  2. Cacluate cause-specific rates by dividing deaths of a given cause into person-years of exposure.
  • This is equivalent to multiplying the overall rate by the ratio of deaths of a given cause to the total.

Steps

  1. Cause-specific probability of death (qi): \(q \cdot dx_i / dx\)
  2. Cacluate cause-specific rates by dividing deaths of a given cause into person-years of exposure.
    • This is equivalent to multiplying the overall rate by the ratio of deaths of a given cause to the total.
lt_ <- 
  lt_ %>% 
    # Cause specific probability of death. 
    mutate(qi = q * dx_i / dx) %>% 
    mutate(Rd = (dx - dx_i) / dx) %>% 
    mutate(md = m * Rd) %>% 
    mutate(mi = m - md) %>% 
    mutate(dxd = dx - dx_i)

age dx dxd dx_i q qi md mi
0 553.221 544.221 9 0.00553 0.00009 0.00546 0.00009
1 38.180 37.180 1 0.00038 0.00001 0.00037 0.00001
2 23.160 22.160 1 0.00023 0.00001 0.00022 0.00001
3 17.689 16.689 1 0.00018 0.00001 0.00017 0.00001
4 14.209 13.209 1 0.00014 0.00001 0.00013 0.00001
5 13.213 12.213 1 0.00013 0.00001 0.00012 0.00001

Modeling CVD-Deleted Mortality

  • We can also apply the methods covered earlier to the cause-deleted life table data!
ages_     <- lt_$age[lt_$age<100 & lt_$age>=25]
deaths_   <- lt_$dxd[lt_$age<100 & lt_$age>=25] 
exposure_  <- lt_$lx[lt_$age<100 & lt_$age>=25]

mort_fit_CVDdeleted <- MortalityLaw(
    x  = ages_,
    Dx  = deaths_,   # vector with death counts
    Ex  = exposure_, # vector containing exposures
    law = "HP2",
    opt.method = "LF2")

Basic CVD Model

G Healthy Healthy Healthy->Healthy CVD CVD Healthy->CVD r_HS Dead Dead Healthy->Dead r_HD(t) CVD->CVD Dead(CVD) Dead (CVD) CVD->Dead(CVD)    r_DCVD(t)         CVD->Dead  r_HD(t) Dead(CVD)->Dead(CVD) Dead->Dead

Basic CVD Model

  • Healthy Death based on modeled (Heligman-Pollard) cause-deleted mortality.
G Healthy Healthy Healthy->Healthy CVD CVD Healthy->CVD r_HS Dead Dead Healthy->Dead r_HD(t) CVD->CVD Dead(CVD) Dead (CVD) CVD->Dead(CVD)    r_DCVD(t)         CVD->Dead  r_HD(t) Dead(CVD)->Dead(CVD) Dead->Dead

Basic CVD Model

  • CVD Death based on age- and cause-specific death rates from cause-deleted life table.
G Healthy Healthy Healthy->Healthy CVD CVD Healthy->CVD r_HS Dead Dead Healthy->Dead r_HD(t) CVD->CVD Dead(CVD) Dead (CVD) CVD->Dead(CVD)    r_DCVD(t)         CVD->Dead  r_HD(t) Dead(CVD)->Dead(CVD) Dead->Dead

Basic CVD Model

  • We also allow for a CVD Death transition based on non-CVD causes.
G Healthy Healthy Healthy->Healthy CVD CVD Healthy->CVD r_HS Dead Dead Healthy->Dead r_HD(t) CVD->CVD Dead(CVD) Dead (CVD) CVD->Dead(CVD)    r_DCVD(t)         CVD->Dead  r_HD(t) Dead(CVD)->Dead(CVD) Dead->Dead

Age-Dependent Transition Matrices

Construct a Markov Trace

trace <- 
  sim_cohort(params)
trace[["natural_history"]][1:10,]
  Healthy     CVD CVDDeath    Dead
0  100000     0.0   0.0000    0.00
1   90391  9506.2   0.2945  103.00
2   81702 18088.0   1.1404  208.58
3   73847 25833.6   2.4853  316.88
4   66745 32822.6   4.2817  428.11
5   60324 39126.6   6.8544  542.45
6   54519 44809.4  11.5751  660.14
7   49270 49930.7  17.3944  781.42
8   44526 54544.1  23.8195  906.57
9   40236 58697.0  31.3741 1035.87

Good News

  • Modeled CVD and non-CVD mortality closely matches the cause-deleted life table data!

Your Turn